clear;
syms theta;
y=sin(theta)+2.0/3*sin(theta)*sin(theta/2)*sin(theta/2);
y_max=double(max(subs(y,theta,[0:0.0005:pi])))
Cr=4.0/y_max